/**
   Copyright 2008 Kirill Ignatiev

    This file is part of kerr-wavefronts.

    kerr-wavefronts is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    kerr-wavefronts is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with kerr-wavefronts.  If not, see <http://www.gnu.org/licenses/>.
 */
#include <iostream>
#include <fstream>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#include <cmath>
#include <getopt.h>
#include <cassert>
using namespace std;
#include "settings.h"
#include "front.h"

#define loop(i,N) for (int i = 0; i < N; ++i)

int window;
const int window_title_length = 1024;
char window_title[window_title_length] = "Wavefronts";

double ortho_x = 410, ortho_y = 240;
double true_ortho_x = 410, true_ortho_y = 240;
double seconds_per_frame = 1./30.;
double simulation_time_scale = 4.;
double current_time = 0.;
vector<wavefront> front;
vector<ic_param*> front_ic;
vector<front_history*> front_hist;
bool paused = false;

void reshape(int width, int height) {
  height = max(height, 1);
  glViewport(0, 0, width, height);

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  double m = min(width/ortho_x, height/ortho_y);
  double s = 0;
  double w1 = width/m, h1 = height/m;
  gluOrtho2D(-w1/2-s, w1/2+s, -h1/2-s, h1/2+s);
  true_ortho_y = h1+2*s;
  true_ortho_x = w1+2*s;

  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();

  glClearColor(0.1, 0.1, 0.4, 0);
  glShadeModel(GL_SMOOTH);
  glEnable(GL_POINT_SMOOTH);
  glEnable(GL_LINE_SMOOTH);
  glEnable(GL_BLEND);
  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  glHint(GL_LINE_SMOOTH_HINT, GL_DONT_CARE);
}

void mouse(int button, int state, int x, int y) {
  if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
    double gX,gY,gZ;
    double model[16], proj[16]; GLint viewport[4];
    glGetDoublev(GL_MODELVIEW_MATRIX, model);
    glGetDoublev(GL_PROJECTION_MATRIX, proj);
    glGetIntegerv(GL_VIEWPORT, viewport);
    gluUnProject(x,y,0, model,proj,viewport, &gX,&gY,&gZ);
    printf("Position (%f,%f); radius %f.\n", gX,gY,sqrt(gX*gX+gY*gY));
  }
}

void timer_display(int value) {
  if (!paused) {
    current_time += seconds_per_frame * simulation_time_scale;
    loop(i,front.size()) front[i].evolve(current_time);
    glutPostRedisplay();
  }
  glutTimerFunc((int)(1000*seconds_per_frame), &timer_display, 0);
}

void keyboard(unsigned char k, int x, int y) {
  switch(k) {
  case 'q': case 'Q': exit(0); break;
  case ' ': paused = !paused; break;
  case '-': simulation_time_scale /= 1.05; break;
  case '+': simulation_time_scale *= 1.05; break;
  case 'a': Show_axes = !Show_axes; break;
  case 's':
    ortho_x /= 1.1; ortho_y /= 1.1;
    reshape(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT));
    break;
  case 'S':
    ortho_x *= 1.1; ortho_y *= 1.1;
    reshape(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT));
    break;
  case 'p': Wavefront_precision *= 1.1; break;
  case 'P': Wavefront_precision /= 1.1; break;
  case 'r': case 'R':
    front = vector<wavefront>();
    loop(i,front_ic.size())
      front.push_back(wavefront(front_ic[i], *front_hist[i]));
    current_time = 0;
    paused = false;
    break;
  default: break;
  }
  glutPostRedisplay();
}

void describe_controls() {
  printf("\
Keyboard controls:\n\
  q,Q    Exit\n\
  space  Pause/unpause\n\
  -,+    Decrease/increase simulation speed\n\
  a      Show/hide axes\n\
  s,S    Zoom in/out\n\
  p,P    Increase/decrease wavefront smoothness\n\
  r,R    Reset the simulation to time t=0.\n\
  left click     Print the coordinates of the mouse pointer.\n\
");
}

void display() {
  snprintf(window_title, window_title_length, "Wavefronts (time %.2f / speed %.1f) (range %.0f,%.0f)", current_time, simulation_time_scale, true_ortho_x/2, true_ortho_y/2);
  glutSetWindowTitle(window_title);
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

  if (Show_axes) {
    glLineWidth(1.);
    glBegin(GL_LINES);
    glColor4d(1,1,1,0.5);
    double w = true_ortho_x, h = true_ortho_y;
#define V glVertex2d
    V(-w/2,0); V(w/2,0); V(0,-h/2); V(0,h/2);
    glEnd();
  }
    
  loop(i,front.size()) front[i].render();

  glFlush();
  glutSwapBuffers();
}

// Options

void print_help(char* progname) {
  printf("\
Usage: %s [OPTION]...\n\
Display the time-evolution of a family of null geodesics in a Kerr geometry.\n\
\n\
Geometry arguments:\n\
     --type=TYPE  Use a geometry of type TYPE, which is schw for\n\
                    Schwarzschild geometry, or kerr for Kerr geometry.\n\
     --charge=Q   A black hole parameter (charge).\n\
     --ang=A      A black hole parameter (angular momentum in units of mass).\n\
\n\
Initial set-up arguments:\n\
     --ic=IC        Use the shape of the null geodesic family.\n\
                      Can be a line, or a circle.\n\
     --distance=D   A parameter for the geodesic family\n\
                      (how far away from the centre it starts).\n\
     --size=S       A parameter for the geodesic family (angle for the circle\n\
                      type (in units of 36 degrees),\n\
                      and length for the line type).\n\
\n\
Simulation arguments:\n\
 --history=H\n\
 --delay=D         Each wavefront is drawn at current time and \n\
                     H-1 previous times, spaced as D.\n\
 -n,--number=N     N wavefronts wil be simulated. Individual options for each\n\
                     wavefront should be given in order.\n\
\n\
Miscellaneous arguments:\n\
 --draw-wavefront-lines\n\
 --draw-wavefront-points\n\
 --draw-history        Draw the space between consecutive wavefronts.\n\
 --correct-brightness=(true|false)\n\
                  Scale brightness of line elements by how much the\n\
                      line elements are expanded from the beginning.\n\
 --intensity=I    The scale of the brightness (e.g., 5000)\n\
 --precision      The larger the precision, the smoother the wavefront.\n\
 -h, --help       Print this help message.\n\
 -u, --usage      Print typical options.\n\
\n\
Keyboard commands:\n\
  Type s,S to zoom in and out.\n\
  Type p,P to increase or decrease smoothness.\n\
  Type a to display axes.\n\
", progname);
}

void print_usage(char* progname) {
  printf("\
Usage: %s --type=[kerr|schw] --ic=[line|circle] -n N --history=H --delay=D\n\
",progname);
}

int main(int argc, char** argv) {  
  static struct option long_options[] = {
    {"type", required_argument, 0, 't'},
    {"ic", required_argument, 0, 'i'},
    {"charge", required_argument, 0, 'q'},
    {"ang", required_argument, 0, 'a'},
    {"distance", required_argument, 0, 'd'},
    {"size", required_argument, 0, 's'},
    {"help", no_argument, 0, 'h'},
    {"usage", no_argument, 0, 'u'},
    {"history", required_argument, 0, 8},
    {"delay", required_argument, 0, 9},
    {"number", required_argument, 0, 'n'},
    {"draw-wavefront-lines", optional_argument, 0, 0},
    {"draw-wavefront-points", optional_argument, 0, 0},
    {"correct-brightness", optional_argument, 0, 0},
    {"draw-history", optional_argument, 0, 0},
    {"precision", required_argument, 0, 10},
    {"intensity", required_argument, 0, 11},
    {"false-intensity", optional_argument, 0, 0},
    {0,0,0,0}
  };

  string hole_type="kerr";
  vector<string> ic_type;
  vector<double> distance, size;
  double charge = 0, angmom = 0;
  int history_length = 1;
  double delay = 4;
  int number = 1;

  int option_index = 0;
  int c;
  while(1) {
    c = getopt_long(argc, argv, "hun:", long_options, &option_index);
    if (c == -1) break;
    switch(c) {
    case 't': hole_type = optarg; break;
    case 'i': ic_type.push_back(optarg); break;
    case 'q': charge = atof(optarg); break;
    case 'a': angmom = atof(optarg); break;
    case 'd': distance.push_back(atof(optarg)); break;
    case 's': size.push_back(atof(optarg)); break;
    case 'h': print_help(argv[0]); exit(0); break;
    case 'u': print_usage(argv[0]); exit(0); break;
    case 8: history_length = atoi(optarg); break;
    case 9: delay = atof(optarg); break;
    case 'n': number = atoi(optarg); break;
    case 0:
      if (!optarg || string("true") == optarg) 
        Sets(long_options[option_index].name, 1);
      else if (string("false") == optarg)
        Sets(long_options[option_index].name, 0);
      else
        printf("Do not know what to do with %s=%s (should true/false).\n",
               long_options[option_index].name, optarg);
      break;
    case 10: Wavefront_precision = atof(optarg); break;
    case 11: Wavefront_intensity = atof(optarg); break;
    case '?': break;
    default: break;
    }
  }
  
  if (optind < argc) {
    cout << "Unexpected extra argument: " << argv[optind] << endl;
    exit(1);
  }

  if (ic_type.size() == 0) ic_type.push_back("circle");

  assert(0 < Gets("precision", 1));
  assert(hole_type == "schw" || hole_type == "kerr");
  loop(i,ic_type.size())
    assert(ic_type[i] == "line" || ic_type[i] == "circle");
  assert(0 <= charge && charge <= 1);
  assert(0 <= angmom && angmom <= 1);
  loop(i,distance.size())
    assert(0 <= distance[i]);
  loop(i,size.size())
    assert(0 <= size[i]);
  assert(1 <= history_length);
  assert(0 < delay);

  double Cdistance = 4, Csize = 10;
  string Cic_type = "line";
  loop(i, number) {
    if (i < distance.size()) Cdistance = distance[i];
    if (distance.size()==0) Cdistance = 10.0 * (i+1);
    if (i < size.size()) Csize = size[i];
    if (i < ic_type.size()) Cic_type = ic_type[i];

    front_history* h = new front_history(history_length, delay);
    // insert front
    ic_param* ic;
    if (hole_type == "schw")
      if (Cic_type == "line")
        ic = new line_param(Cdistance,Csize);
      else
        ic = new circle_param(Cdistance,Csize/36*M_PI);
    else {
      kerr_param p(angmom/2, charge); // units of mass
      if (Cic_type == "line")
        ic = new kerr_line_param(Cdistance,Csize,p);
      else
        ic = new kerr_circle_param(Cdistance,Csize*36/180*M_PI,p);
    }
    front_ic.push_back(ic);
    front_hist.push_back(h);
    wavefront w(ic, *h);
    // Colours
    static double colours[7][3] = {
      {1,0.5,0.5}, {0.5,1,0.5}, {0.5,0.5,1},
      {1,1,0.5}, {1,0.5,1}, {0.5,1,1},
      {0.8,0.8,0.8} };
    int j = i%7;
    loop(c,3) w.colourmap[c] = colours[j][c];

    front.push_back(w);
  }

  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
  glutInitWindowSize(450,360);
  window = glutCreateWindow(window_title);
  glutDisplayFunc(&display);
  //glutIdleFunc(&display);
  glutMouseFunc(&mouse);
  glutReshapeFunc(&reshape);
  glutKeyboardFunc(&keyboard);
  glutTimerFunc((int)(1000*seconds_per_frame), &timer_display, 0);

  describe_controls();
  glutMainLoop();
  return 0;
}
